Show the code
pacman::p_load(tidyverse, tmap, sf)Hands-on Exercise 7c: Analytical Mapping
1 About this Exercise
In this exercise, we will learn how to plot analytical maps.
2 Getting Started
2.1 Installing and Loading R Packages
For this exercise, other than tmap, we will use the following packages:
- tidyverse for tidying and wrangling data
- sf for handling geospatial data
The code chunk below uses p_load() of pacman package to check if the abovementioned packages are installed in the computer. If they are, they will be launched in R. Otherwise, pacman will install the relevant packages before launching them.
2.2 Importing the Data
For the purpose of this hands-on exercise, a prepared data set called NGA_wp.rds will be used. The data set is a polygon feature data.frame providing information on water point of Nigeria at the LGA level.
We will use read_rds() function to import the data into R.
Show the code
nga <- read_rds("data/rds/NGA_wp.rds")
glimpse(nga)Rows: 774
Columns: 9
$ ADM2_EN <chr> "Aba North", "Aba South", "Abadam", "Abaji", "Abak", …
$ ADM2_PCODE <chr> "NG001001", "NG001002", "NG008001", "NG015001", "NG00…
$ ADM1_EN <chr> "Abia", "Abia", "Borno", "Federal Capital Territory",…
$ ADM1_PCODE <chr> "NG001", "NG001", "NG008", "NG015", "NG003", "NG011",…
$ geometry <MULTIPOLYGON [m]> MULTIPOLYGON (((548795.5 11..., MULTIPOL…
$ total_wp <int> 17, 71, 0, 57, 48, 233, 34, 119, 152, 66, 39, 135, 63…
$ wp_functional <int> 7, 29, 0, 23, 23, 82, 16, 72, 79, 18, 25, 54, 28, 55,…
$ wp_nonfunctional <int> 9, 35, 0, 34, 25, 42, 15, 33, 62, 26, 13, 73, 35, 36,…
$ wp_unknown <int> 1, 7, 0, 0, 0, 109, 3, 14, 11, 22, 1, 8, 0, 37, 88, 1…3 Basic Choropleth Maps
We will plot
total <- tm_shape(nga) +
tm_polygons("total_wp",
palette = "Greens",
style = "equal",
lwd = 0.1,
alpha = 1) +
tm_layout(main.title = "Distribution of All Water Points",
main.title.size = 1,
legend.outside = FALSE)
total
functional <- tm_shape(nga) +
tm_polygons("wp_functional",
palette = "Blues",
style = "equal",
lwd = 0.1,
alpha = 1) +
tm_layout(main.title = "Distribution of Functional Water Points",
main.title.size = 1,
legend.outside = FALSE)
functional
tmap_arrange(total, functional, nrow = 1)
4 Choropleth Map for Rates
4.1 Deriving Proportion of Functional Water Points and Non-Functional Water Points
nga <- nga %>%
mutate(pct_functional = round(wp_functional/total_wp,2)) %>%
mutate(pct_nonfunctional = round(wp_nonfunctional/total_wp,2))
glimpse(nga)Rows: 774
Columns: 11
$ ADM2_EN <chr> "Aba North", "Aba South", "Abadam", "Abaji", "Abak",…
$ ADM2_PCODE <chr> "NG001001", "NG001002", "NG008001", "NG015001", "NG0…
$ ADM1_EN <chr> "Abia", "Abia", "Borno", "Federal Capital Territory"…
$ ADM1_PCODE <chr> "NG001", "NG001", "NG008", "NG015", "NG003", "NG011"…
$ geometry <MULTIPOLYGON [m]> MULTIPOLYGON (((548795.5 11..., MULTIPO…
$ total_wp <int> 17, 71, 0, 57, 48, 233, 34, 119, 152, 66, 39, 135, 6…
$ wp_functional <int> 7, 29, 0, 23, 23, 82, 16, 72, 79, 18, 25, 54, 28, 55…
$ wp_nonfunctional <int> 9, 35, 0, 34, 25, 42, 15, 33, 62, 26, 13, 73, 35, 36…
$ wp_unknown <int> 1, 7, 0, 0, 0, 109, 3, 14, 11, 22, 1, 8, 0, 37, 88, …
$ pct_functional <dbl> 0.41, 0.41, NaN, 0.40, 0.48, 0.35, 0.47, 0.61, 0.52,…
$ pct_nonfunctional <dbl> 0.53, 0.49, NaN, 0.60, 0.52, 0.18, 0.44, 0.28, 0.41,…functional_rate <- tm_shape(nga) +
tm_polygons("pct_functional",
n = 10,
palette = "Blues",
style = "equal",
lwd = 0.1,
alpha = 1,
legend.hist = TRUE) +
tm_layout(main.title = "Rate Map of Functional Water Points by LGAs",
main.title.size = 1,
legend.outside = TRUE)
functional_rate
5 Extreme Value Maps
5.1 Percentile Map
5.1.1 Data Preparation
nga <- nga %>%
drop_na()percent <- c(0, 0.01, 0.1, 0.5, 0.9, 0.99, 1)
var <- nga["pct_functional"] %>%
st_set_geometry(NULL)
quantile(var[,1], percent) 0% 1% 10% 50% 90% 99% 100%
0.00 0.00 0.22 0.48 0.86 1.00 1.00 5.2 Creating get.var function
get.var <- function(vname, df) {
v <- df[vname] %>%
st_set_geometry(NULL)
v <- unname(v[,1])
return(v)
}5.3 Percentile Mapping Function
percentmap <- function(vnam, df, legtitle=NA, mtitle="Percentile Map"){
percent <- c(0,.01,.1,.5,.9,.99,1)
var <- get.var(vnam, df)
bperc <- quantile(var, percent)
tm_shape(df) +
tm_polygons() +
tm_shape(df) +
tm_fill(vnam,
title=legtitle,
breaks=bperc,
palette="Blues",
labels=c("< 1%", "1% - 10%", "10% - 50%", "50% - 90%", "90% - 99%", "> 99%")) +
tm_borders() +
tm_layout(main.title = mtitle,
title.position = c("right","bottom"))
}percentmap("total_wp", nga)
6 Box Map
ggplot(data = nga,
aes(x = "",
y = wp_nonfunctional)) +
geom_boxplot()
boxbreaks <- function(v,mult=1.5) {
qv <- unname(quantile(v))
iqr <- qv[4] - qv[2]
upfence <- qv[4] + mult * iqr
lofence <- qv[2] - mult * iqr
# initialize break points vector
bb <- vector(mode="numeric",length=7)
# logic for lower and upper fences
if (lofence < qv[1]) { # no lower outliers
bb[1] <- lofence
bb[2] <- floor(qv[1])
} else {
bb[2] <- lofence
bb[1] <- qv[1]
}
if (upfence > qv[5]) { # no upper outliers
bb[7] <- upfence
bb[6] <- ceiling(qv[5])
} else {
bb[6] <- upfence
bb[7] <- qv[5]
}
bb[3:5] <- qv[2:4]
return(bb)
}get.var <- function(vname,df) {
v <- df[vname] %>% st_set_geometry(NULL)
v <- unname(v[,1])
return(v)
}var <- get.var("wp_nonfunctional", nga)
boxbreaks(var)[1] -56.5 0.0 14.0 34.0 61.0 131.5 278.0boxmap <- function(vnam, df,
legtitle=NA,
mtitle="Box Map",
mult=1.5){
var <- get.var(vnam,df)
bb <- boxbreaks(var)
tm_shape(df) +
tm_polygons() +
tm_shape(df) +
tm_fill(vnam,title=legtitle,
breaks=bb,
palette="Blues",
labels = c("lower outlier",
"< 25%",
"25% - 50%",
"50% - 75%",
"> 75%",
"upper outlier")) +
tm_borders() +
tm_layout(main.title = mtitle,
title.position = c("left",
"top"))
}tmap_mode("view")
boxmap("wp_nonfunctional", nga)